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Abstract 

In this paper, we highlight the topological properties of leader neurons whose existence 
is an experimental fact. 

Several experimental studies show the existence of leader neurons in population bursts 
of 2D living neural networks [1,2]. A leader neuron is, basically, a neuron which fires at 
the beginning of a burst (respectively network spike) more often that we expect by looking 
at its whole mean neural activity. This means that leader neurons have some burst triggering 
power beyond a simple statistical effect. In this study, we characterize these leader neuron 
properties. This naturally leads us to simulate neural 2D networks. To build our simulations, 
we choose the leaky integrate and fire (IIF) neuron model [3,4], which allows fast simulations 
[5,6]. 

Our IIF model has got stable leader neurons in the burst population that we simulate. 
These leader neurons are excitatory neurons and have a low membrane potential firing 
threshold. Except for these two first properties, the conditions required for a neuron to be a 
leader neuron are difficult to identify and seem to depend on several parameters involved in 
the simulations themself. However, a detailed linear analysis shows a trend of the properties 
required for a neuron to be a leader neuron. Our main finding is: A leader neuron sends a 
signal to many excitatory neurons as well as to a few inhibitory neurons and a leader neuron 
receives only a few signals from other excitatory neurons. 

Our linear analysis exhibits five essential properties for leader neurons with relative im- 
portance. This means that considering a given neural network with a fixed mean number 
of connections per neuron, our analysis gives us a way of predicting which neuron can be 
a good leader neuron and which cannot. Our prediction formula gives us a good statistical 
prediction even if, considering a single given neuron, the success rate does not reach hundred 
percent. 
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1 Introduction 

Disassociated in vitro rat brain neuron cultures show, under a variety of experimental contexts, 
a spontaneous electrical activity. This activity manifests itself by a rapid succession of igni- 
tions of a large fraction of neurons named collective bursts (or network spikes) [7-10]. This 
bursting activity can be measured, for example, by multi-electrode array methods [1,2, 1 1] or by 
fluorescence [12]. The study of the initiators of these bursts is one of the conceptual problems 
underlying the spontaneous electrical activity. A recent study [1] showed that some "first to fire" 
cells exist. Later [2], through a detailed analysis of data obtained by the multi-electrode array 
methods, established that some cells are triggering bursts beyond a simple statistical effect of 
being the first. These particular cells are called leaders. The analysis [2] indicates that the long 
term dynamics of the leaders is relatively robust, evolving with a half-life of about one day. [2] 
also shows that these leaders are not only the main initiators of bursts, but, the burst itself carries 
traces (or hints) indicating which of the leaders has initiated the burst. In this respect, one can 
view the culture as an amplifier of the signal emitted by the leaders. 

It is clear that the scenario described above calls for a theoretical explanation. One of the 
ways pursued is the one of (bootstrap) percolation in random networks [13] which gives insights 
into the spread of the initial ignition and distinguishes between localized and de-localized spread- 
ing. In this picture, the burst itself is viewed as the "giant component" of the percolation process, 
and this picture is verified from many different angles [12]. 

In this paper, we address the delicate question of what makes a neuron a leader. Is it stim- 
ulation? Activity threshold? Special connectivity properties of the network? The natural way 
of reaching this aim theoretically is to simulate a random neural network. Our results show that 
leadership is an effect of a combination of several ingredients which can be quantified by a simple 
relation between several natural parameters of the neuron. We obtain these results by a simple 
model of a randomly connected network with the usual leaky integrate and fire [3] mechanism 
for ignition of neurons. This well known model [4] is tested in numerical event-based simula- 
tions that provides, in principle, unbiased simulations [14]. This model only takes into account 
simple properties of neurons (integrate and fire) [5,6] but it suffices to reproduce basically all the 
findings of [2]. Thus we trade realism for conceptual simplicity, but in this kind of investigations 
the complexity comes from the neural network itself and not from the neuron model. 

The results of our simulations can be presented as follows: to each neuron one can assign 
a leadership score and neurons with the highest scores are leaders. It turns out that leaders can 
be characterized as being excitatory neurons and having a low membrane potential firing thresh- 
old. Apart from these two first properties, the conditions required for a neuron to be a leader 
are difficult to identify and depend on several parameters involved in the simulations themself. 
Indeed, a detailed linear analysis shows a trend of the properties required for a neuron to be a 
leader. The main finding of this paper is a formula for the leadership score that exhibits five es- 
sential properties for leaders with relative importance. This formula gives a very good statistical 
prediction. Therefore, we can conclude that leadership is a random effect, and that leaders are 
formed naturally from a balanced combination of inputs, outputs, their local neighborhood and 
their own properties. Basically a leader sends a signal to a lot of excitatory neurons as well as to 
a few inhibitory neurons and a leader receives only a few signals from other excitatory neurons. 
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2 The simulations 

In this section we explain how we construct our leaky integrate and fire simulations and which 
kind of parameters we use. 

2.1 Building the simulations 

A single neuron is already a complex biological object. But, basically, neurons are electrically 
excitable cells that are composed of a soma (cell body), a dendritic tree and an axon. Dendrites 
and axons connect to each other (through synapses) to create a neural network. Thus, to make 
a model, we can, for example, choose a geometry and some rules that produce the rule of con- 
nection between neurons [15] or construct a random neural network with a fix mean number of 
connections per neuron. 

Our simulations are built as follows. First we construct a matrix of synaptic weights W that 
represents a random neural network of G IN neurons where is a parameter and \/N G IN 

as so we can consider that the neurons are placed on a grid y) G ]N^| x,y < \/N | C 

with periodic boundary conditions. In the Euclidean geometry of R^, we consider that every 
neuron n (n = 1, 2, . . . , A^) is a circle of radius r G R+ (equivalent to the spatial extension 
of the dendrites). Each axon has a length ^ G R+, approximately Gaussian distributed (with 

probability density function close to /(£, L) = ^^^ e vvhere L G R+ is the mean axon 

length and A = ^min — 1, — is the standard deviation (and 1 < L < |v^)). In 
our simulations, the neuron spatial parameters are the dendrite size r and the mean axon length 
L. To every axon, we also associate a direction (respectively an angle) 9 G [0, 27r) uniformly 
distributed (with probability density function g(9) = ). Figure 1 illustrates the network 
geometry. 

To construct the matrix W G MAr(R), we proceed as follows: If the axon (with finite length 
i, without width and angle 9 (respectively direction)) of the neuron m {m = 1, . . . , N) intersects 
the circle (dendrites) of the neuron n ^ m {n = 1, . . . , N) then the connection exists (i.e.Wnm 7^ 
0), otherwise Wnm = 0. The biological knowledge tells us that the connection between neurons 
is oriented (electrical impulses move from axons to dendrites). Thus we have a non-symmetric 
matrix of synaptic weights W representing the neural network. 

A fix parameter w G R+ named the synaptic weight is the value we give to Wnm when the 
connection exists {i.e.Wnm 7^ 0). Note that, in reality, some neurons are inhibitors. For these, 
we use a negative synaptic weight — w G R_. In our simulations the proportion r G [0, 1] of 
inhibitory neurons is also a parameter. 

To simplify the terminology, if Wnm 7^ 0, we will say that the neuron n is the son of the 
neuron m and the neuron m is the father of the neuron n. One given neuron can have more than 

'Negative axon lengths are not allowed and we impose a maximum length of for each axon. With this 

last condition, no axon can cover all the network in any direction. The rationale for this choice is that, in the 
experiments [2], probably no axon covers all the test tube in a given direction. 
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one father and more than one son. Of course, a given neuron can be both, father and son of 
another neuron. 



cb cb cb 



e 



Figure 1 : A 3 x 4 neurons rectangle. Each neuron is a 
circle with a radius r G 1R+ (dendrites), in this figure 
r < 1. The axon of one of the neurons is drawn. 
We see that this neuron is connected to the neuron 
on the left of the bottom row. Note that the axon 
direction is given by ^ € [0, 27r). Note also that if 
r > 1, unlike in this figure, then we automatically 
have nearest neighbour connections. 



Once the matrix of synaptic weights is constructed we assign two real values to each neuron 
n: V* the membrane potential firing threshold and Vn{t) the membrane potential which depends 
on the time t E R+. Without loss of generality, we can consider that the membrane potential ^ 
oscillates roughly between and {V*) = 1 (where (V*) is the mean of all the neuron membrane 
potential firing thresholds). The time evolution of the membrane potential V"„(t) is computed 
using the leaky integrate and fire [3] neuron model. 

The dynamics is the following, every neuron n (respectively every membrane potential Vn{t)) 
which integrates the electrical signal it receives from the other neurons is able to fire (when 
Vn{t) > V*) and then sends a signal to the other neurons. Just after firing the membrane potential 
is reset to 0. The membrane potential also decreases with a leaky time r (loss of memory). 

More precisely, the membrane potential evolution equation of the neuron n is the evolution 
equation of a leaky integrate and fire neuron model [3] 



N 



Vnit + dt) = Vn{t)e-- (1 

Loss of memory 



Reset after firing 



m=l 



Vr„, t 



dr. 




Connections 



where 



dt is the infinitesimal time step, 

r E is the leaky time (or characteristic time of the membrane [3]), 

Xn[Ki(^)] is, for the neuron n, the indicator function of the set of potentials bigger or equal 
to V*. Namely, Xn[Vn{t)] = 1 whenever Vn{t) > V* and Xn[Vn{t)] = otherwise. The 
distribution of all the characteristic membrane potential firing threshold is approximately 

Gaussian with a probability density function close to h (V*) 
AV*, the standard deviation, is a parameter. 



1 



AVW21T 



(V*-\ 
e 2AV* 



where 



^We know that the membrane potential oscillates more or less between —70 and —50 mV [3]. But to simplify, 
we will approximately scale the membrane potential between and 1. 
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• dnm = dmn IS the distance between the neuron n and the neuron m. In the experiments 
related in [2], one finds a 2D density of about 2000 This leads us here to fix the metric 
scale in our simulations: If we observe an expected mean number of neurons of 2000 
neurons inside the grid square, we assume that this square measures 1 mm^. Nevertheless, 
we can also measure the distance in an arbitrary unit given by the expected number of 
neurons in a row but we need the metric scale to compare our signal propagation velocity 
to the experimental one, 

• f is the propagation velocity of the signal in the neural network, 

• Bt^n is, for each time t and each neuron n, a Gaussian random variable with mean a > 
and standard deviation |. Roughly speaking Bt^n ~ A/" (cr, |) for all n and for all t > 0. 
To generate this Gaussian noise, we use a Box-MuUer transform, 

• ^It^n is a Markovian exponential clock with a mean rate u. For most t, one has fi^ .„ = 0, 
but sometimes l^t „ = 1- More precisely, the distribution of the time intervals At of the 
continuous set t where l^t „ = follows the exponential distribution ^e~^ {to is the mean 
value). Finally, it means that Bt n-^t,n is a noise „ which acts on the membrane potential 
of neuron n only when the exponential clock fi^ „ rings. 

By looking at (1), we see that, in principle, we can compute every membrane potential all the 
time. But it is not useful and too complicated to compute every membrane potential continuously. 
Choosing a fixed time step 5t > also introduces some problems like synchronized events and 
systematic errors [16]. In the program, made in Perl, we compute the membrane potential only 
when an event happens^. At this time, we check if the neuron fires or not. So, in the program, we 
store the events to find out the next membrane potential to update. In other words, our simulations 
(1) are event-based simulations. The firing time of neurons is not discretized but computed event 
per event at the machine precision level. This way of doing provides, in principle, unbiased 
simulations [14]. 

Note that if "every" neuron fires at (almost) the same time (this may happen in a burst) then 
all the membrane potentials are reset to at the same time; the neurons are synchronized. To 
avoid that this situation occurs too often, when we initialize a sequence, we give to each neuron n 
a membrane potential in [0, V*) (uniformly distributed). Note that we reinitiate a sequence each 
time a burst is too long, so it is possible that during a simulation the neurons get synchronized. 

Note also that if the mean exponential clock rate uj is too small relatively to the leaky time 
r, it does not really matter which kind of distribution we use for the noise i?t „. This means that 
if w ^ r then we can use Bt^n = for all t and for all n as the noise distribution B^^n and the 
result of the simulations remains the same. 

^Here, by event we mean the moment when a neuron receives a signal from another neuron or when a neuron 
receives some noise (/.e.the Markovian exponential clock Q,t,n rings). 
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2.2 Values of parameters and simulation choices 

The parameters are fixed in a way to observe activity like the one observed in the experiments [2] 
when we restrict the observation to sixty neurons. This section provide a list of the order of 
magnitude of the parameters used in (1). Note that the reader can find almost all the information 
in [17]. 

• We made simulations for different sizes G IN. But beyond ~ 4 it does not affect the 
presence of bursts and leaders. This means that we already get a leader when we consider 
only four neurons. 

• In fact fixing the dentrites size r, the mean axon length L and the axons length distribution 
fix the distribution function for the sons (out-degree) and the mean number of connections 
per neuron (noted (connections/neuron)). Remark that the number of sons per neuron 
distribution looks like the axons length distribution while the number of fathers per neuron 
distribution is always approximately Gaussian distributed (data not shown). Note also that 
the geometry we choosed for our simulations induces that the successful pre-burst shows 
locality, like the one observed by the authors in [2]. Finally in most of our simulations, we 
chose r = 0.85 and L = 7, then (connections/neuron) = 11. 

• The synaptic weight w is surely a function of the neuron and of the time (the neuron's 
learning). As soon as we choose mean number of connections per neuron and considering 
the kind of neural activity we want, we can fix w. We use \n j^. Note that if w = (V*) = 
1 (where (V*) is the mean of all the neuron membrane potential firing thresholds) then all 
neurons are too strongly correlated and if w = then all neurons are independent. 

• For the inhibitory proportion r, the choice of r = 0.2 is reasonable. Depending on the 
neural culture type, the inhibitory cells proportion can change approximately from 20% to 
30% (see [12]). 

• In our simulations, the leaky time r is about 100 ms. Note that the order of magnitude is 
the correct one [18]. To be exact we compute with r = 0.1 s, this means that the used time 
unit is the second. 

• The propagation speed v in a neural network is between 50 and 100 see [17]. In our 
case, V should be replaced by the propagation speed in the axon which is bigger. However, 
we have decided to use this range throughout. 

• The standard deviation of the membrane potential firing threshold AV* is also a parameter. 
We used < AV^* < 0.2{V*), which guarantees that neurons are quit similar and still the 
membrane potentials do not get synchronized too often. 

• For the exponential clock rate we use cu between 0.001 and 0.01 s. 

• The mean noise a is determined by the kind of neural activity we want to have (observe). 
If there is too much noise then we loose the locality in the successful pre-burst because we 
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have a lot of uncorrelated activities. Note that the mean membrane potential (without con- 
nection) V = limT-s>oo ^ Jq V (f) di is about"^ We want the mean membrane potential 

V to be high enough to get some neural activity due to the noise. But we also want that 

V + v\i < (V*), so that one spike does not necessarily create a burst. 

One might think that to keep the neural activity one can vary the values tu and a keeping 
the ratio - constant because the mean membrane potential V is about — . But to preserve 
the neural activity, the mean membrane potential V is not that relevant. However it is 
necessary to keep constant the mean time the membrane potential of neurons need to reach 
the value V*. From this point of view the relation between a, cu and w is a first passage 
time problem. Note also that if the value of the mean membrane potential V is too low 
then there is no neural activity due to the noise (this means no neural activity at all). To 
conclude, even if the exact value of the mean membrane potential V is not fundamental, 
this value must not be too far from the mean of all neurons membrane potential firing 
threshold {V*) (see [3]). 

The parameters N, r, L, AV* and r are the neural network parameters and the parameters r, v, 
w, a and cu are the dynamical parameters. 





Neural network parameters 


Dynamical parameters 


Known 


r = 20% (inhibitory proportion) 
r, L, (mean number of connections) 


r ^ 100 ms (leaky time) 
V ^ 100 Sfi (propagation speed) 


To be set 


N (computing-memory) 
AV* (firing thi^eshold standard deviation) 


(V*) — V > \N > (synaptic weight) 
cr, LU (noise parameters) 



Table 1: Known and unknown parameters of the simulations. The unknown parameters will be set in 
Section 4.1. 



3 General appearance of the simulations and some definitions 

In this section we show what is the product of our simulations and we give the definitions of what 
we call leader neuron, according to [2]. 

3.1 General appearance 

Like in the experiments [1,2, 10], we want to look at the neural activity. The product of our 
simulations consists of a list of ordered spikes times. For each spike i, the list gives us U, the 
time of the i^^ spike, and Ui, the neuron that fired. Then, considering only the spike chronology, 
we analyze the "spike sequence" of a subset (or all) neurons of our simulations. 

'*The mean membrane potential charge due to the noise is ct + ae^^ + ae^'^^ + . . . = — ^ 12: jf iii jg 
small. 
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15 [ms]: 



Isolated Aborted pre-burst Trigger Successful pre-burst Burst 




Figure 2: Picture of the temporal evolution of a simulation for the whole neural activity with the param- 
eters: N = 900, r = 0.85, L = 7, r = 0.2, Ay* = 0.05, = 5 ^, w = 0.1, r = 100 ms, a = 0.008, 
uo = 0.001 s, nburst = 500, (5tburst = 15 ms and ^tisoiated = 5 ms. The color code, that gives a repre- 
sentation of the different definitions explained in the following section, is the following: Isolated spike, 
aborted pre-burst, trigger, successful pre-burst, and burst. 

3.2 Definitions 

Figure 2 represents a part of a simulation and leads us to define sequences in the simulations. 
This section is roughly the same than [2, Section 2.3]. Precisely we will define mathematically 
what we call a burst and a leader and we will also explain what we call mutual information. In 
order to make the present paper comparable to [2], we adhere strictly to the same definitions. So 
the reader who red [2] or knows these definitions can skip this part and move directly to the next 
section. 

3.2.1 Definition of bursts and their triggers 

First we divide all the spikes of our list into four classes: burst, successful pre-burst, aborted pre- 
burst, isolated spike, and each spike belonging to one and only one class. The precise definitions 
are given below, but basically a spike is in a burst if it is in a group of many spikes that follow 
each other closely. If a spike is not in a burst it can be in a successful pre-burst or aborted 
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pre-burst. But it must be anyway is in a sequence of spikes that are close enough in time so that 
communication between them is still possible. The distinction between successful or aborted pre- 
burst depends on whether the spike sequence is eventually followed by a burst or not. Finally, all 
other spikes are isolated. 

More precisely, we use three parameters, ribmstj bursts and (5tisoiatcd- We first look for in- 
terspike gaps of length at least (5tisoiated, and we divide the set of all spikes into disjoint subsets 
TZk of consecutive spikes, where A; is a running index. The TZk have the property that if (the 
spike) i and i + 1 are in 7?.fc then U+i — U < 5tisoiated, while if i E TZk and i + 1 G 7?.fc+i then 
ti+i—ti > (5tisoiatcd- The rationale is that if the interspike gap is so large that the spike z + 1 is due 
to noise and not to the axon signaling of the precedent (in time) spike i (i.e. tj+i — U > 5tisoiatod 
then the spikes belong to separate subsets). 

We now want to further divide each TZ^ into the burst itself, characterized by a high density 
of spikes, and its precursor which immediately precedes it in time but has a lower density. Each 
of the TZk contains at least one spike, but may contain many. In each TZk we first look for a spike 
that is followed by at least riburst — 1 spikes in TZk within a lapse of time 5t burst- Denoting this 
spike's index by = i^{k), this is the first index in TZk with the property: 

^i*+nburst-l ~ ^i* ^ '^^burst • (2) 

If the condition (2) is met then we subdivide TZk into two disjoint sets at i*: TZk = 'Pk^ CPk 
can be empty). The index up to the last index in TZk make up the burst Bk. 

If the condition (2) is never met in TZk, then the set TZk is not subdivided. This is called an 
aborted pre-burst if it has more than one spike and is an isolated spike otherwise. 

For these k where such an can be found, we check if it is also the first spike in TZk- If so, 
this burst is an immediate burst that has no successful pre-burst, so that TZk = ^k- For all the 
other bursts we divide TZk = 'Pk^ ^k where T^k are the index i G TZk with i < i^,, and Bk are the 
others. The letter P refers to successful pre-bursts and B refers to bursts. 

Between bursts the neural activity is much lower, and we label these periods as quiet. Each 
quiet period, denoted Qt, for a running index h, is actually a concatenation of consecutive TZk^ 
that did not contain a burst in the earlier stage, or else an empty set. We now renumber the set of 
index as follows 

{1, . . . , M} = Qi U Pi U i3i U Q2 U P2 U i32 U . . . , 

with Bi, the fe'th burst, the corresponding successful pre-burst (if it exists else is empty) 
and Qb the corresponding quiet phase. We define the trigger of burst h as the neuron rii^, where 
ib is the number of the first spike in Pj, (resp. Bi, if is empty like for an immediate burst). 



N 


Analyzed neurons 


'^burst 


5tburst [ms] 


'^^isolatcd 




60 


20 




2 


900 


15 


900 


> 100 


3 V 



Table 2: Burst detection parameters. 

The parameters riburst, ^^burst, and 5tisoiated that we used for = 900 are shown in Table 2. 
Note that these parameters can easily be adapted for other values of the number of neurons A^. 
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Remark also that the time parameters 5tisoiatcd and 5tburst have to be smaller than the character- 
istic time of the membrane r. The rationale for these choices are: Firstly, to cross the distance 
of the biggest possible axon '^VN, the signal, at the velocity v, needs a time of This com- 
putation gives us the order of (5tisoiatcd- This means that only the singnaling velocity is taken into 
account to compute (5tisoiatcd- Secondly, even in a burst, the neural activity is not regular, so to 
balance this effect we need (5tburst to be high enough. Finally, for nburst we want it high enough 
to guarantee that a burst includes many neurons. 



3.2.2 Leaders 

For every successful pre-burst, aborted pre-burst or immediate burst, the neuron which fires first 
is called the trigger. Since we are interested in their special role, we first need to make sure that 
triggers are not just neurons with high neural activity, which are statistically more often the first 
ones to fire. The following discussion will exhibit that leaders are over-proportionally more often 
triggers. 

To qualify a neuron as a leader, we require that a trigger's probability to lead bursts should 
be significantly higher than its probability to fire. Let M be the total number of bursts in a 
simulation. For each neuron n, we define a„ as the number of times that neuron n has spiked 
during the simulation. Now, let us consider a spike, and term by Qn the probability^ that neuron 
n has fired that spike: g„ = On/ ctn'- The probability for the neuron n to be a spurious, or 
random trigger F times in M bursts is given approximately^ by the binomial distribution Pn{F): 



M 



Q^ii-quf'-^. (3) 



In the limit of large M and reasonable g„ this distribution is approximated by a Gaussian of mean 
Mqn and variance Mg„ (l — Qn)- On the other hand, we denote by /„ the actual number of bursts 
that neuron n leads (note that J2n fn = 

Thus we have a scale on which to test triggering. We define a„, a "leadership score", and 
decide that a neuron is a leader if it scores at least 3 standard deviations above the natural expec- 
tation value. The criterion for leadership of neuron n is therefore 

fn-Mqn 

an = — , > 3 . (4) 

Because we expect a ~ A/'(0, 1) to be Gaussian, (4) corresponds to a p-value of 0.001. 



3.2.3 Mutual Information 

We estimate the mutual information of two neurons (pair of neurons) in a series of time intervals 
based on empirical probabilities, according to [19]. Using the division (of time) into aborted 

^unbiased estimator of the considered probability 

^It would be exactly the binomial distribution if all the spikes were independent. All the spikes are certainly not 
independent because neurons are connected. So the binomial distribution is not the correct one but gives a good 
enough approximation (see [15]). 
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pre-burst, successful pre-burst and burst intervals, we define a firing event in a given interval as 
h = 0,1, where stands for no spike and 1 stands for at least one spike. We denote by Nn{h) the 
number of times neuron n had an event h'ma series of A^int intervals^. The probability assigned 
to the event h is thus^ Pn{h) = Nn{h) /N-mt- The number of joint events in which neuron n had 
an event hi and neuron n' had an event h2 is given by Nn,n'{hi, h2). The joint probability of 
events hi, /12 for neurons n, n' is then given by^ Pn,n'(^i, ^2) = Nn,n'{hi, h2)/Ni^t. The mutual 
information between two neurons is then given by 



Note that the mutual information between two neurons is if the two neurons fire independently 
and is In 2 if the two neurons fire always together (see [20]). These two extreme values give us 
the variation scale of the mutual information for a pair of neurons. 



In this section we explain our results. We prove that our simulations have stable leaders and 
we identify the topological properties of these leaders. An analysis of the first to fire neurons 
exhibits that the leader properties we find are the correct ones. 

4.1 About the parameters and some statistics 

The next step consists in analyzing our simulations knowing that our results will be more con- 
sistent if we span as many parameters as possible. However, our goal is not to obtain a transition 
phase diagram [21], so we are not going to vary all parameters (most of the them were fixed in 
Section 2.2 through Table 1). Up to now we stayed unclear about some parameters, including 
dynamical parameters like w, a and uj. So we still have to find (and fix) a range for the parame- 
ters AV^*, w, cr and a; in a way to obtain "reasonable" simulations^ which involves distinguishing 
bursts, successful pre-bursts and quiet periods. One of the best ways to do that mathematically 
is to look at the time interspike distribution and compare with the experimental one (data not 
shown, see [15]). The analysis of the membrane potential profile and the time length distribution 
are also possible in the simulations (see [15]). 

4.2 Existence of leaders 

In Section 3.2.2, we defined, for every neuron n a leadership score (see (4)) that compares the 
triggering frequency to the leadership expected frequency (compute with respect to the neural 
activity frequency). Figures 3 show the leadership score as a function of the time for a few 

^Note that iVjnt = M for the burst intervals and TVint < M for the successful pre-burst intervals. 
^Here by "reasonable" simulations, we mean to observe a neural activity more or less Uke the one in Figure 2 
but without a precise predicting spike timing [22, 23]. 




?ii,/i2e{o,i} 




Pn,n'{hi, /la) 
Pn{hl)Pn'{h2) ' 



4 Results 
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neurons. In these figures we see that our simulations have leaders. Note that this result is very 
robust. We mean that the choice of the used parameters in simulations is not important. As soon 
as we use "reasonable"parameters, our simulations have leaders. This is still true when we only 
look at a few neurons instead of all them (Figure 3 a). Even more, we can change the axons length 
distribution and the simulations still have stable leaders. Note also that even a network with a 
uniform probability of connections^ shows bursts and leaders but, in this case, we do not observe 
any spatial locality. 




(a) The leadership score when we look at only 60 neu- 
rons Uke in [2] . We see that we have leaders and even 
very good ones. We also see that the average lead- 
ership score {a) is about 0. Note that in this graph, 
each point represents 600 seconds of simulation that 
contains more than 2500 bursts. Remark that the best 
leader among 60 neurons is not necessarily a leader 
while we consider all the neurons (see Figure 3b). In 
the restricted case, the neurons detected as leaders can 
be first to fire neurons (see Section 4.4). 



(b) The leadership score when we consider all the 900 
neurons. We can say that these neurons are the "true" 
best leaders of the simulation. Note also that, luckily, 
one of the best leaders was in the 60 selected neurons. 
But even when it is not the case, when we consider 
only 60 neurons, some neurons show very good lead- 
ership scores. On this graph, each point represents 
about 2000 seconds of simulation that contains more 
than 8800 bursts. 



Figure 3: Leadership score of a few neurons for the simulation shown in Figures 2. Note that the leadership 
score depends on the number of neurons considered for the analysis (see neuron 634, 643 and Section 4.4). 

In our simulations, with > 400 about 5% to 10% of the neurons have a mean leadership 
score (4) bigger than the threshold 3 when we consider all the neurons. While considering only 
60 neurons, about 15% neurons have a mean leadership score bigger than the threshold 3. Also 
note that even if the leadership score can fluctuate during the simulations, its mean value stays 
stable (z.e.its standard deviation is small). Of course, the leader ratio decreases when the mean 
noise a increases. This fact is logical because in a too noisy neural activity it becomes more 
difficult to identify the trigger of each burst. Note also that it seems that when AV* (standard 
deviation of the membrane potential firing threshold) is between 0.05 and 0.1 the leader ratio 
seems maximal (data not shown). 



'^This means a network that is not build with the kind of geometry showed in Figure 1 . 
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In the experiments [2], one electrode captures in general the spike signal of several neurons 
[2, Figure la]. But, in our simulations, up to this point we considered each neuron separately, 
as if one electrode measures only one neuron. Remark that the leadership score of a group of 
neurons also reach the value 3 for some groups [15]. 

At this point, our simulations reproduce basically all the findings of [2] (presence of leaders, 
signature of bursts, force of leaders,. . . ) see [15]. The advantage of simulations is that we have 
access to a lot of details that we did not have access to in the experiments and we can analyze the 
topology of the network to find out the properties of leaders. 

4.3 Leader properties 

4.3.1 Some facts about the leader properties 



Firing threshold compared to leadership score 




"^■^ 10 20 30 40 50 

Leadership seore (ct-n) 

(a) The membrane potential firing threshold {V*) com- 
pared to the leadership score (an). We see that good 
leaders have a low membrane potential firing threshold 
{{V*) = 1). Parameters: N = 900, r = 0.85, L = 7, 
r = 0.2, AF* = 0.1, v = h-^,\N = 0.1, t = 100 ms, 

ms 

a = 0.008, uj = 0.001 s, nturst = 100, Jiburst = 15 
ms and St^soUtcd = 5 ms. 



Triggering frequency compared to leadership score 




10 20 30 

Leadership score (a^) 



(b) The triggering frequency (/„/M) compared to the 
leadership score (q!„). These two things are propor- 
tional. Parameters: N = 900, r = 0.85, L = 7, 
r = 0.2, AV* = 0.05, v = 5 ^, w = 0.1, 
T = 100 ms, cr = 0.008, u ^ 0.001 s,"nburst = 100, 
<5iburst = 15 ms and (5iisoiated = 5 ms. 



Figure 4: Some properties of leaders for simulations with parameters close to the one we showed in Figure 
2. Note that each point in these figures represents an average of the results for 5 simulations using the 
same neural network and the same parameters. Each simulation contains more than 10000 bursts and the 
simulated time was longer than 25 minutes. 

In Figures 4 we see clearly that the leaders seem to have some special properties compared 
to the other neurons. More precisely. Figure 4a shows that good leaders typically have a low 
membrane potential firing threshold. This fact exhibits that to be a good leader, a neuron must 
have the property of firing early. But, this property is not the only one a neuron needs to have 
to be a leader. Indeed the best leader in Figure 4a is not the neuron with the lowest membrane 
potential firing threshold. Thus, another property of the leaders is easy to highlight: the leaders 
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are excitatory neurons^" (data not shown). In Figure 4b, we see that, like in [2], the better a leader 
is, the more bursts it triggers. 

We can also have a look at the leadership score compared to the neural activity ratio of the 
neurons shown in Figure 5. Firstly, remember that to compute the leadership score we use the 
neural activity ratio q (see (4)). Secondly, look at Figure 5 and remark that there are as much 
leaders with a high neural activity ratio as leaders with a low neural activity ratio. So, clearly 
and as expected, for a given neuron there is no link between having a high leadership score and 
the neural activity ratio, even if a good leader triggers a lot of bursts (see Figure 4b). 

Before trying to find more leader properties, we add more remarks about Figure 5. We note 
that there exist some neurons with a very high positive leadership score (a > 10) but there is no 
neuron with a very low leadership score {a < — 10). We also note that if the standard deviation of 
the membrane potential firing threshold /SV* is then most of the neurons have a neural activity 
ratio near q = {m Figure 5g=g^ = 0.11%as expected). 



Leadership score compared to the neural activity ratio 



0.001 0.002 0.003 

Neural activity ratio {q-n) 



Figure 5: The leadership score {a) compared to the neu- 
ral activity ratio {q). Note that, as before, each point 
in this figure represents an average of the results for 5 
simulations using the same neural network. Parameters: 
N = 900, r = 0.85, L = 7, r = 0.2, AV* = 0.1, u = 5 
5^, r = 100 ms, w = 0.1, a = 0.008, uj = 0.001 s, 
^burst = 100, 5tburst = 15 ms and dtisoiaXed = 5 ms. 



Unfortunately, the other relations between a property and the leadership score that we tried 
to highlight did not exhibit an as good correlation as the relations showed in Figures 4. We also 
note another important fact: In the special case where we use a network with every membrane 
potential firing threshold fixed to {V*) {i.e. /SV* = 0) or/and without inhibitory neurons (i.e. r = 
0), our simulations still have stable leaders. These two facts lead us to the following hypothesis: 
the properties of leaders are a non trivial combination of properties. 

To find these properties, we perform the following experiment: We take a given network^ ^ 
and make several simulations with this network but change some of the dynamical parameters (r, 

'"Remark that in the special case where w ^ t, we rarely observed a few inhibitory neurons with a leadership 
score bigger than 3. This can happen only when these inhibitory neurons have a very low membrane potential firing 
threshold, a few excitatory fathers and a lot of inhibitory fathers. The reason is as follow, because a; <C r, the 
network gets synchronized after the first burst and these particular inhibitors have the faculty to fire even before the 
true trigger of the burst. Beside the type of their fathers implies a very low neural activity, so their leadership scores 
are potentially good. Note that the few top best leaders are always excitatory neurons. 

"This means that we choose N, r, L, p, AV*, r and build a realization of a network {i.e. W and are known 
for all n). 
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V, w, a and u) in a way of keeping a reasonable simulation. Then, we compute the leadership 
score of each neuron and we do this for each simulation. As a result, we find that the best leader 
is not always the same even if most of the time a leader stays leader. Figure 6 illustrates this 
fact and this forces us to admit that the leader properties depend on the dynamical parameters. 
Nevertheless, we can probably find a trend of the properties of a typical leader by doing a linear 
least squares analysis for all our simulations. Indeed, in Figure 6 more than 90% of the leaders 
keep on being leaders and non leaders keep on being non leaders under the change of dynamical 
parameters. 

Figure 6: The leadership score of two different simu- 
lations using the same neural network. The red dashed 
lines separate leaders and non-leaders for both simula- 
tions. We observe that the best leader is not necessarily 
the same and some neurons are leaders with a set of pa- 
rameters but are not leaders with the other set. We also 
observe that most of the time a good leader continues 
being a good leader in both simulations. Like for Fig- 
ures 4, the error bars are computed by averaging over 5 
simulations using the same network and dynamical pa- 
rameters w, a and w. Common parameters: = 900, 
r = 0.85, L = 7, r = 0.2, AF* = 0.1, v = b 
j^, r = 100 ms, nburst = 100, 5tburst = 15 ms and 
<^iisoiated = 5 ms. 

Before going through this linear least squares analysis, note that in the experiments related 
in [2], the leaders change during the days spent in vitro with a typical life time of about one day 
(see [2, Figure 3a]). While, in the simulations, the leaders stay the same in a given network if we 
do not change the dynamical parameters of the simulation. So, one of the interesting facts in the 
simulations is that to change the leaders we do not need to change the network, some changes 
in the dynamical parameters are enough (we consider the synaptic weight w as a dynamical 
parameter). 

4.3.2 Finding the leader properties 

To find out the leader properties we do the following linear least square analysis. First, we build 
a vector a = {an}n=i n '^hat contains the leadership score a„ of each neuron n of a given 
simulation. Then we construct a matrix P = {Pni}i<:n<N i<j<6 ^^^^ contains all the neuron 
properties (excitatory or inhibitory neuron and membrane potential firing threshold) and all the 

'^Here by a reasonable simulation, we mean that we want to observe quiet periods, successful pre-bursts and 
bursts in a way to find leaders. Note that this range of reasonable parameters is quite large and the relation between 
the dynamical parameters is non trivial, see remark in Section 2.2 and Section 4. 1 . 



Leadership score comparison 
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Leadership score a„(w = 0.1, ct = 0.0081, u = 0.001 ! 
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connection properties at the first order . More precisely: P„i = 1 if the neuron n is an excitatory 
neuron and —1 otherwise. P„2 = ^"av* ^ Al^* > and otherwise. P„3 is number of 
excitatory neurons connected to the axon of neuron n (this means the number of excitatory sons 
of the neuron n, see Section 2.1) and P„4 the number of inhibitory sons. Identically, is the 
number of excitatory fathers of the neuron n and the number of inhibitory fathers. 

Our linear least square method, like in [24], consists in finding a vector x = {xi}-^^ g such 
that the euclidean norm 

||Px — a|| (5) 

is minimized. This particular x gives us the relations (at the linear order) between the first order 
properties'"* of a neuron and its leadership score. 

By repeating the calculation (5) for each of our simulations and averaging over all our simu- 
lations, we obtain the typical solution X = (1.22, —2.33, 0.13, —0.19, —0.13, 0) . This means 
that a good prediction p„ for the leadership score a„ of the neuron n in a network is given by 

Pn = 1.22 ■ P„i - 2.33 ■ P„2 + 0.13 ■ P„3 - 0.19 ■ P„4 - 0.13 ■ P„5 , (6a) 
1.22 ■ excitatory,, - 2.33 ■ ^""^r^ + 0.13 ■ # excitatory sons„ 

ifA\/*>0: Pn= AV J n ^^^^ 

—0.19 • # inhibitory sons^ — 0.13 ■ # excitatory fathers^ 

where (excitatory^ =) P„i = 1 if the neuron n is an excitatory neuron and —1 otherwise, 
Pn2 = ^"av* ^ AV* > and otherwise, P„3 (= # excitatory sons„) is the number of 
excitatory sons (of neuron n), Pni (= # inhibitory sons,„) is the number of inhibitory sons and 
PnX) (= # excitatory fathers^J is the number of excitatory fathers. 

Note also that (6) do not depend on the axons length distribution. Even more, by building a 
fully random network (without axons and dendrites, i.e.no particular geometry) and fixing only 
the mean number of connections per neuron, we obtain the same formula. 

Precisely, equations (6) mean that a typical leader in our simulations is an excitatory neuron, 
that it has a low membrane potential firing threshold and that it has a lot of excitatory sons but a 
few inhibitory sons and a few excitatory fathers. Equations (6) in this exact form are valid for a 
mean number of connections per neuron of about eleven ((# connections/neuron) = 11). 

Now, considering a given neural network (that we did not use to obtain (6)) and using (6), 
we can predict, even before running our simulation which neuron will be a leader and which 
not. Figure 7a shows a typical prediction. Remark that this prediction is pretty reliable because: 
All the neurons with a prediction j9„ lower than are not leaders and all the neurons with a 
prediction pn bigger than 6 are leaders > 3). More precisely. Figure 7b shows the prediction 
efficiency'^ Pe as a function of the leadership score a. In Figure 7b we can also remark that 
the prediction is correct for more than 90% of the neurons (pe(a) > 0.9 Va). This means that 



'^We call connection properties at the first order all the direct connections of the neuron in the network contained 
in the matrix of synaptic weights W (this means the number of sons and fathers (see Section 2.1)). 

'"'We call first order properties all the neuron properties and all the connection properties at the first order 
^^We call prediction efficiency the coefficient Pe{a) ~ — ^'i^ng ™"'"^'p'|"'"^"°"" y^jjere correct prediction^ is 1 if 
the neuron n is a leader a„ > 3 (respectively not a leader < 3) in the simulation and its prediction p„ > 3 
(respectively p„ < 3), otherwise: correct prediction„ = 0. 
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Figure 7b shows that the trend of prediction is correct. Of course Figure 7b also shows that there 
are some mistakes in the predictions especially for prediction (p„) between and 6 (see Figure 
7a). It is important to remember that the leaders are function of the dynamical parameters and 
our coefficients in (6) (which are only computed at the linear order) are not. 



Leadership score compared to prediction 



Prediction elBciency as a function of the leadership score 




^ 0.7 



Prediction (pji ) 



Leadersiiip score (a) 



(a) The leadership score (a„) compared to the predic- 
tion (pn) for a particular simulation. The red dashed 
lines separate leaders and non-leaders for both the sim- 
ulation and the prediction. As usual, the error bars of 
each point in this figure are computed using 5 simula- 
tions done with the same network and parameters. 



Figure 7: Effectiveness of the prediction formula (6). Parameters: = 900, r = 0.85, L = 7, r = 0.2, 
AV* = 0.1, u = 5 ^, r = 100 ms, w = 0.1, a = 0.0081, u = 0.001 s, riburst = 100, <5tburst = 15 ms 
and (Jiisoiated = 5 ms. 

Figure 8 shows the prediction efficiency^^ as a function of the leadership score a. Most 
of the simulations presented in Figure 8 show that the trend of prediction of (6) is correct 
(Pe(a) > 0.9, Va). Precisely, for a range of simulations with different standard deviation of 
the membrane potential firing threshold (AV*) (Figure 8 simulations B, D, E and F), for dif- 
ferent mean noise (a) and exponential clock rate {to) (Figure 8 simulations C and F) the trend 
of prediction is correct. So by looking only at the simulations B, C, D, E and F of Figure 8, 
we can say that (6) give a pretty good prediction (even without adapting the coefficient of (6) 
with the dynamical parameters). But for the simulation A the trend of prediction is less good 
(Pe(tt) < 0.9, Va > 0). The fact is the following, for the simulation A, the inhibitory propor- 
tion r = 0. This means that for simulation A we did not stay in the parameters range that we 
gave in Section 2.2 and that we used to find out (6). We can say that the fact that (6) does 
not give a very good trend of predictions for the simulation A is a kind of edge effect. With 
Xr=o = (-0.8927, -3.23, 0.07, Ai, -0.02, As)^ in (6) (where Ai, As G R), we obtain a much 
better prediction for simulation A (pe(«) > 0.92, Va). Remark that the first number (—0.8927) 



(b) The prediction efficiency of Figure 7a. For 
each neuron n we attach a number, named correct 
prediction,!, which is 1 if the neuron n is a leader 
q;„ > 3 (respectively not a leader a„ < 3) in the simu- 
lation and its prediction p„ > 3 (respectively p„ < 3), 
otherwise: correct prediction„ = 0. For each value of 
the leadership score a, we compute the prediction effi- 

/ \ ^ \ ^ conect prediction , . , . 

ciency Pe[a) = ' j which is plot 

as a function of the leadership score a. 
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(the coefficient of P„i = excitatory^) of Xr=Q has no physical meaning because all neurons are 
excitatory neurons (if r = then P„i = excitatory^ = 1 \fn). For the same reason Ai, A2 in 
Xr=o can be any real number because P„4 (=# inhibitory son„) and P^g (=# inhibitory fathers^) 
are always if all neurons are excitatory neurons. All this means that the relative importance of 
the leader properties change by changing the parameters of the simulation but the trend in (6) 
is correct. The sign of each coefficient in (6) is robust (as soon as the coefficient has a physical 
meaning). 



Prediction efficiency as a function of the leadership score 















— A 






— B 






— C 






— D 






— K 






— F 











Leadership score (a) 



Figure 



8: The prediction efficiency 

correct prediction 



Pe{a) = 

^ as a function of the leader- 
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ship score a for different simulations (see Figure 7b). 
A) r = 0, w = 0.0666, AV* = 0.1, a = 0.008 and 

io = 0.001 s, B) r = 0.2, w = 0.1, AV* = 0.05, 
cr = 0.008 and uj = 0.001 s, C) r = 0.2, w = 0.1, 
AV* = 0.1, a = 0.014 and lo = 0.002 s, D) r = 0.2, 
w = 0.1, AV* = 0, a = 0.008 and u = 0.001 s, 
E) r = 0.2, w = 0.1, AV* = 0.2, a = 0.007 and 
u = 0.001 s, F) r = 0.2, w = 0.1, Ay* = 0.1, 
(7 = 0.008 and uj = 0.001 s. Common parameters: 
N = 900, r = 0.85, L = 7, t; = 5 r = 100 



ms, nburst 
ms. 



100, 5tburst = 15 rns and 5t 



isolated 



We can also explain why we do not use a bigger matrix P. Considering a given neural 
network, we can extract a lot of different properties for all neurons n: the number of fathers, 
the number of sons, the number of fathers 'fathers (second order property'^), the number of loop 
order 2, 3, 4, 5 . . . in which the neuron n participates,. . .We can also reduce the network to an 
excitatory network and so on. We can as well try to improve our results by using a lot of tricks 
and/or cut off. But using the properties at the first order seem to be an easy and logical way of 
doing. Especially considering two facts: 1) A son will not necessarily fire after one of his fathers 
did, so what to say about the sons of this son? 2) The leadership score depends on the dynamical 
parameters. Taking that fact into account improves more the result than considering the second 
order properties in the network. 

To conclude, (6) cannot be taken as a general law and (6) is only a good approximation of the 
leadership score in the range of parameters we chose in Section 2.2. However, we pretend that 
the signs in (6) are robust under the change of parameters. This means that our typical leader 
profile is the right one. 

We pretend that even if (6) are computed on small scale simulations, the trend of the typical 
leader profile is the correct one. Even more we propose to use the following leader profile if the 



'^This means that we need to compute to extract this property. 
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mean number of connections per neuron vary 



ifr,AV*>0: p„ 



1.22 ■ excitatory^ - 2.33 ■ 

42 * inhibitory sons^^ 

r (# connections /neuron) 



AV* 

1.14- 



1.14 



# excitatory sons^ 



( 1 — r ) (# connections / neuron) 



# excitatory fatliers 
(1— r)(# connections /neuron) 



(7) 



This equation was computed knowing that in (6): (# connections/neuron) = 11. 

We have already explained in Section 4.2 that we observe leaders in our simulations as soon 
as we use reasonable parameters. Further more, considering a given simulation, we can store 
the neurons using the leadership score. After that, if we decide to remove some of the (best) 
leaders from the network and then run again our simulation we still observe leaders (as good as 
before). Finally, in a leaky integrate and fire neuron model there are leaders in the networks and 
the neurons which fit better with the typical leader profile are the leaders. 

4.4 First to fire neurons 

In this section we characterize the first to fire neurons. From this point, we will call first to fire 
neurons the neurons that fire in a short time after the trigger. This means that the first to fire 
neurons fire roughly in the successful pre-burst. 



# fathers/sons as a function of tlie position in the burst 
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Figure 9: The mean number of fathers/sons per neuron 
as a function of the time position of a spike of this neu- 
ron in the burst. The time bin in which we average the 
number of fathers/sons per neuron is 1 ms. The time po- 
sition is reserved for the trigger of the burst. The mean 
number of connections per neuron in the network (noted 
(# connections /neuron)) is about 11. The typical num- 
ber of fathers of an active neuron /„ is about 11.5 (see 
(8) for definition). We obtained this figure by averaging 
over more than 35000 bursts. Parameters: N = 900, 
r = 0.85, L = 7, r = 0.2, AV* = 0.05, v = 5 
5^, w = 0.1, r = 100 ms, a = 0.008, uj = 0.001 s, 
"burst = 150, Jtburst = 15 ms and (Jtisoiatcd = 5 ms. 



Figure 9 shows the mean number of fathers (or sons) per neuron as a function of the time 
position of a spike of this neuron in the burst. To obtain Figure 9 we average over all the bursts 
of a given simulation without taking into account if the trigger of the burst was a leader or not. 
Nevertheless, when we look at the mean number of sons per trigger'^ we remark that the trigger 
of the burst has, on average, more sons than their followers. This fact is consistent with our 
previous analysis. Remember that Figure 4b shows that leaders trigger more bursts than other 
neurons. This means that the trigger is, on average, a leader. Remember also that (6) tells that a 
leader has typically many excitatory sons (this means a lot of sons as well). 

''The time position in Figure 9 is reserved for the trigger of the burst. 
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Figure 9 also shows that the first to fire neurons have, on average, more sons than their 
followers and have, on average, less sons than the trigger. This means that very often several 
leaders are present in the successful pre-burst (this is linked with the triggering force of leaders, 
data not shown but see [15]). So the fact that first to fire neurons have typically a lot of sons is 
consistent with (6). 

In Figure 9 the mean number of sons per neuron tends to the mean number of connections 
per neuron in the network. This means that, except fot the trigger and the first to fire neurons, the 
mean number of sons per neuron for the neurons which fire during the burst itself is not relevant. 

Figure 9 shows explicitly that the mean number of fathers per neuron of the first to fire 
neurons is higher than the mean number of connections per neuron in the network. This means 
that the first to fire neurons are characterized by a high number of fathers. This is consistent with 
the fact that neurons with a lot of fathers light up early in the burst, respectively in the successful 
pre-burst (see [25]). In facts neurons with a lot fathers have an higher neural activity than neurons 
with a few fathers (data not shown). And again, this is consistent with [25] and Figure 9 which 
shows that active neurons in the burst are, on average, neurons with a lot of fathers (typically 
more fathers than the mean number of connections per neuron in the network). 

In Figure 9 the mean number of fathers per neuron (unlike the mean number of sons) does not 
tend to the mean number of connections per neuron in the network but seems to increase again 
after the decay of the first to fire neurons. This is consistent with the following fact: in a burst, 
after a while, neurons are able to fire again (as soon as their membrane potentials are charged 
enough). And, of course, the neurons with a lot of fathers are able to fire again sooner than the 
others. 

What about the mean number of fathers of the trigger? Figure 9 shows a tendency of the 
trigger to have more fathers than the mean number of connections per neurons. This last fact 
seems to contradict our predictions (6). But in Figure 9 the mean number of fathers of the trigger 
is approximately what we call the typical number of fathers of an active neuron /„: 



where # fathers^ is the number of fathers of neuron n and g„ is the neural activity ratio of neuron 
n. Figure 9 shows that (except for the first to fire) the mean number of fathers per neuron is 
about fa the typical number of fathers of an active neuron and shows that the number of fathers 
of the typical trigger is not different from fa. These two facts can explain the term concerning 
the number of fathers in (6), but still, we can think it contradicts (6). 

To clarify this conflict, let us have a look at Figure 10. This figure shows, as Figure 9, the 
mean number of fathers (and sons) per neuron as a function of the time position of a spike of 
this neuron in the burst. In Figure 10 the standard deviation of the membrane potential firing 
threshold AV* is much higher than the one in Figure 9. In Figure 10, we see that the number of 
sons of the trigger as well as the one of the first to fire neurons is higher than the mean number 
of connections per neuron in the network. And we see that the number of fathers of the trigger 
is lower than the mean number of connections per neuron in the network. But the number of 
fathers of the first to fire neurons is higher than the mean number of connections per neuron in 



N 




(8) 
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Figure 10: The mean number of fathers/sons per neuron 
as a function of the time position of a spike of this neu- 
ron in the burst. The time bin in which we average the 
number of fathers/sons per neuron is 1 ms. The time 
position is reserved for the trigger of the burst. The 
mean number of connections per neuron in the network 
(noted (# connections /neuron)) is about 11. The typ- 
ical number of fathers of an active neuron f,, is about 
11.9 (see (8) for definition). We obtained this figure 
by averaging over more than 35000 bursts. Parameters: 
N = 900, r = 0.85, L = 7,r = 0.2, AV* = 0.2, v = 5 
5^, w = 0.1, r = 100 ms, a = 0.007, u = 0.001 s, 
"burst = 150, Jtburst = 15 ms and ^tisoiatcd = 5 ms. 



the network. This force us to conclude that the simulation shown in Figure 10 does not contradict 
(6). The reason for that is that the membrane potential firing threshold plays a more important 
role in this case. 

We can also remark that in Figure 9 like in Figure 10, the mean number of fathers per neuron 
tends approximately to the typical number of fathers of an active neuron fa while we look forward 
in time in the burst. 

Now, let us clarify the observation of Figures 3 concerning the neurons 634 and 643. Neuron 
643 is one of the best leaders of the simulation of Figure 3b while restricting the observation 
to 60 neurons only implies that neuron 634, which is not a leader (see Figure 3b), appears as a 
leader (see Figure 3a). The fact is that both neurons 634 and 643 have the profile of first to fire 
neurons (mean number of fathers close to the typical number of fathers of an active neuron). But 
the membrane potential firing threshold of the two neurons explain why only neuron 643 can be 
a leader: V^^^ = 1.009 and V^^^ = 0.869 ((V^*) = 1 and AV = 0.05). So while restricting the 
observation to 60 neurons only, it is possible that the first to fire neurons appear as leaders. In the 
experiments [2] it was impossible to distinguish between true leaders and first to fire neurons. 

Finally we have a competition between two facts: To be a leader, a neuron needs to be able to 
fire early (in order to do that it must have a low membrane potential firing threshold and a lot of 
fathers) and to be a leader, a trigger must trigger more bursts than we expect. Typically a leader 
has a good triggering power (this means a lot of excitatory sons) and a relatively low neural 
activity ratio (this means less fathers than the typical number of fathers of an active neuron). 
That is the message of (6) (respectively (7)). Remember that (6) was obtained by averaging over 
simulations, so (6) express only a tendency of what kind of neurons the leaders are. 



5 Conclusion 



Experimental studies show the existence of leader neurons in population bursts of 2D living 
neural networks [1,2]. This leads naturally to a question: Do leaders also exist in neural network 
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simulations? In this study, we have proved that stable leaders exist in leaky integrate and fire 
neural network simulations. 

In our simulations, we saw that leaders depend on the network but also on the dynamical 
parameters. Because the leaders mainly depend on the network and on the own properties of 
the neurons, we were able to find some important properties for these leaders. Firstly, they are 
excitatory neurons and have a low membrane potential firing threshold (ability of firing early). 
Secondly, they send a signal to a lot of excitatory neurons and to a few inhibitory neurons (trig- 
gering power). Finally, they receive only a few signal from other excitatory neurons (trigger 
more bursts than we expect by looking at the neural activity). 

The study of more biological plausible models [5, 6] (like Hodgkin-Huxley) would be very 
interesting too and could probably also produce stable leaders. These models could probably be 
used to find out other leader properties. 
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